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Abstract. A general method to construct recombinant tree approximations for sto- 
chastic volatility models is developed and applied to the Heston model for stock price 
dynamics. In this application, the resulting approximation is a four tuple Markov pro- 
cess. The first two components are related to the stock and volatility processes and take 
values in a two dimensional Binomial tree. The other two components of the Markov 
process are the increments of random walks with simple values in { — 1, +1}. The result- 
ing efficient option pricing equations are numerically implemented for general American 
and European options including the standard put and calls, barrier, lookback and Asian 
type pay-offs. The weak and extended weak convergence are also proved. 



1. Introduction 

Contrary to many mathematical models, the discrete counterpart of the celebrated Black 
& Scholes [1] model came after its continuous version and it is generally accepted that this 
simple Binomial approximation by Cox, Rox, and Rubinstein [5J has been instrumental in 
the better understanding and the applicability of the model. Rubinstein [25J states that 
"the Black & Scholes model is widely viewed as one of the most successful in the social 
sciences and perhaps, including its binomial extension, the most widely used formula, with 
embedded probabilities, in human history" . 

This widespread use and practicality is extended by further research. In particular, 
stochastic volatility models have been introduced to address the volatility smiles observed 
in option markets and the heavy tails and high peaks of the underlying asset distributions. 
Hull and White [19], Chesney and Scott [5], Stein and Stein [29], Heston [17], and Hagan 
et al. 16 among many others assume a bivariate diffusion framework in which a separate 
stochastic process represents the dynamics of asset price volatility. In all these models, the 
asset price process St and its volatility factor process Y t satisfy the following stochastic 
differential equations 

dS t = S t [fidt + f(Y t )dWt], 
dY t = fi Y \Y t )dt + a Y \Y t )dZ u 

with correlated Brownian motions W, Z and different choices for the functions fi Y (y) , a Y (y) 
and f{y). 

In this paper, we construct an approach that provides a recombining tree approximation 
for all stochastic volatility models of the above type. This approximation as the Cox, Ross 
& Rubinstein (CRR) model easily constructs a discrete time financial market that itself is 
arbitrage free and as such allows for simple analysis of related complex instruments. 
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For specificity, we implement our methodology on the Heston model. As well known, 
among stochastic volatility models it deserves special attention because of its ability to 
provide closed form solutions for European options through Fourier transform. This unique 
feature allows for an efficient and quick calibration of the model to European options. 
However, for derivative products with early exercise features such as American options, 
closed-form solutions do not exist even under the Heston model. Hence, numerical methods 
such as binomial tree, finite difference schemes or Monte Carlo simulation have to be used 
to evaluate American and other exotic options under stochastic volatility models. 

In any market with a non-constant volatility, the CRR methodology encounters a basic 
difficulty. Indeed, since the volatility changes at each time, the nodes do not recombine 
in the lattice and this fact results in an exponential and thus a computationally explosive 
tree that cannot be used in many realistic situations. Nelson & Ramaswamy [24] was the 
first to construct a computationally simple binomial process which approximates a diffusion 
process given in the form 

dY t = n(Y, t)dt + a(Y, t)dZ t . 

They solve the node recombination problem by transforming the process given in the above 
equation into a process X(Y,t) such that the instantaneous volatility of the transformed 
process is constant. Hilliard and Schwartz [18] follows this method to develop binomial 
trees for continuous-time risk-neutralized diffusion processes of a special form. 

Our main tool is to apply correlated random walks in order to approximate diffusion 
processes. A correlated random walk is a generalized random walk in the sense that the 
increments are not identically and independently distributed, but they only satisfy some 
Markov type of conditions. The exact definition will be given in Section [3] These processes 
naturally lie on a grid, and their Markov structure allows for a efficient computation of 
option prices. 

The idea to use correlated random walks for approximating diffusion processes goes back 
to Gruber and Scwhweizer [15] and to Kusuoka [22] . In [15] , the authors prove a convergence 
result for one dimensional diffusion process that satisfy stronger regularity conditions that 
appear in stochastic volatility models. In [22], Kusuoka uses (also in one dimension) an 
original technique to modify random walks in order to get a diffusion in the limit. Again 
the regularity conditions that he assumes are stronger that arising in stochastic volatility 
models. 

Our approach is also similar to of Kusuoka and modify the correlated random walks on 
a multi-dimensional Binomial tree by adding a predictable process times \fh where h is the 
size of the time step. We then use this freedom to choose the predicable process together 
with an appropriate choice of the conditional probabilities to construct a Markov process 
that weakly converge to the stochastic volatility model. This construction is explained in 
Section [3] The weak convergence of our approximation is given in Section [4] Then, the 
approximating martingale measures are constructed so that the modified tree under these 
measures asymptotically match the first two conditional moments. These fact allows for a 
straightforward convergence proof. We also note that this approach was successively used 
by the last two authors [5] to prove convergence of a market with trading costs. 

Our extensive numerical experimentation is reported in our final section. We leave the 
detailed description of the computational studies to that section and here simply state that 
our algorithm works efficiently compared to all existing methods for the Heston model. 

We emphasize that our tool can also be applied for a general type of stochastic volatility 
models and for GARCH models. Indeed, there is extensive research connecting GARCH to 
the stochastic volatility models. We refer the reader to Duan [10] [11] [TSJ [13] , Nelson [23] 
Ritchken & Trevor [27] and the references therein. 

Clearly, there are several other successful computational approaches to stochastic models, 
including the ones based on partial differential equations, semi-analytic methods and Monte- 
Carlo simulations. Here we do not survey all these results but compare our numerical results 
with the appropriate ones in the section that outlines our numerical experimentations. 
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In the literature, tree based methods have also been considered. Beliaeva & Nawalkha 
[5] is the most recent of these studies and we refer the reader to [2] and the references 
therein for these studies. However, our approach differs from these earlier studies in two 
fundamental points. Firstly, our approximation is recombinant by construction while in the 
previous studies recombination is achieved through truncation. Also, our tree is arbitrage 
free and we provide a proof of convergence. 

2. The Heston Model 

Consider the Heston model, 

dS t = St(rdt + ^ t dW t ), 

dv t = n(6 — v t )dt + r)\/v t dWt, 

with initial conditions So,vo > 0, given positive parameters r,K,6,r) and two Brownian 
motions W, W with a constant correlation p € (—1, 1). The constant r > is the interest 
rate and S is the stock price process. As it is standard, we also assume that 

2n6 > j] 2 . 

Then, the Heston equation has a unique positive solution in R^_, see for instance [7]. 

The main goal of this paper is to construct a discrete approximation of this model. For 
this purpose, it is more convenient to work with a transformed system of affine equations 
driven by independent Brownian motions. Therefore, we set 

x t := In St, y t := — - px t , 
V 

so that 

(2.1) dx t = fi x (xt,yt)dt + a(x t ,yi)dWt, 



where 



dyt = n v {x u yt)dt + \/r){l - p 2 ) a(x t ,y t )dB t , 



(J, x (x,y) :=r- ^-rj{y + px), p y (x,y) := — - pr + \ (prj - 2k) (y + px), 
2 Tj 2 

W t - pW t 



B t ~— F=r> a ( x ,y) '■= V(y + P x ) + ' 



and z + = max(0, z). One may directly verify that B is also a standard Brownian motion 
independent of W . 

3. Derivation of the approximation 
We fix a time horizon, or equivalently a maturity, T > and a time discretization 

h := -, 

n 

with a large integer n. We then use two dimensional correlated random walks to approximate 
the diffusion processes given by (|2.ip . Indeed, consider the random walks {A^™\ Y^"' ) } fc=0 
of the form 

k 

(3.1) 4 n) := xo + v^E^' 

i=l 

k 



(3-2) lf } := yo + v^l-P 2 ) E^> 

i=i 

where xq :— ln(so), yo ■— (^o/v) ~ P x o and (C X jC r )' s are random variables with values in 
{ — 1,1}. In the sequel, we always use the initial data 
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We construct a probabilistic structure so that the four tuple (XjP ,Y^ n forms 
a Markov chain weakly approximating the solution of (12.1[) . To achieve this we also need 
to introduce a modification of this discrete Markov chain. Indeed, for given predictable 
processes a, j3, we introduce 

(3.3) X| n) := xW + y/hj&kg, 

(3.4) Y k (n) := F fe (n) + - P 2 ) 4 fc = l,...,n. 
Clearly, the convergence of (X, F) is equivalent to the that of (X, F) as 

H^-W- jr^H = 0(Va), ||y(")_yW|| =0(v^), 

where for any exponent k, we use the standard notation 0(h k ) to denote a generic random 
variable of the order h k and o(h k ) denotes a random variable that converges to zero after 
divided by h k . 

Our goal is to construct a sequence of probability measure P( n ) and stochastic processes 
q,(™),/3(«) such that 

{(X[nt/T]> Y [nt/T])}t=a i( X t' 2/t))S=0> 

where denotes weak convergence. We provide the definitions in the next section. 

In view of the martingale convergence Theorem 7.4.1 in P3], to establish this convergence, 
it is essentially sufficient to match the first and the second conditional moments. Indeed, 
for a positive integer k, set 

and let e£T l '[ ■ ] be the conditional expectation • [Tk] with respect to the probability 

measure p("). Then, the moment matching conditions are the following equations, 

(3.5) Ejft[*£ n) " = M^-i^-i)/* + °W 
(3-6) W-^[]=^W 1I l£l)fc + o(/i) 
(3-7) Ej& [(.*<"> - l^) 2 ! = ^(I^iljft + o(*) 



(3-8) [(Y^ - ^-i)"J = ^(1 - P>Wi, ^-i)^ + 

We also need conditions on the the covariances. However, since (W 7 , B) in (|2.1[) arc 
independent, this condition is simply reduced to the requirement that £^ and are con- 
ditionally independent given J-'k—i- 

Observe that we need to solve four equations and the number of unknowns or parameters 
to choose are four as well; the corrections a,/3 and two probabilities, 

(3-9) Pk := f£\(& = 1), q k := Pjft(# = !)• 

This construction would provide a financial market which is asymptotically arbitrage 
free. However, a slight modification of the above procedure would also ensure that each 
discrete market itself is free of arbitrage. In our model, the discrete stochastic process 

{exp(-rkh) exp(X ( k n) )} k=0 , 

is the approximation of the discounted price process. Hence, we replace the first order 
condition (I3.5|) by requiring that above process is a martingale, i.e., 

(3-10) Et ) 1 [cxp(-r/ l )exp(4" ) )-exp(^ ) 1 )] =0. 

In fact, (|3.5[) and (|3.10l) are asymptotically equivalent and both would be sufficient to prove 
convergence. However, in our numerical experimentation we observe that this modification 
substantially improves the accuracy of the resulting numerical algorithms. 
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We continue by constructing p(") and cS n \ fi^ 71 ' satisfying the equations (|3.10[) . and 
[3.8)1 . Indeed, by (|3.10p we directly calculate that 

(1 + a k )^\[^} - & k -i$_i = o(h). 

Hence, 

(1 + a^)^-!)^-!^]^-! = K-i) 2 + o(h). 
We use this and calculate that 

E& ((*<"> - X^f) = r,h ((1 + a k f - (4"_\) 2 + o(h)) + o(h). 
We assume that a k is varying slowly in fc, which we later verify. Hence, 

((^i n) - 4" ) 1 ) 2 ) = ^(1 + 2d fc ) + o(h). 
We now compare the above equation with (|3.7I) to conclude that 

l + 2a^a 2 (^"J 1 ,y fe ^ 1 ) + (/ l ). 
Using (|3.6[) and (|3.8p . we obtain the same equation for /3. Hence, we conclude that 

afe = Pfc = 2 h 

We use the above identity and the freedom on the order o(h) to define the processes a,f3 
below. The below definition contains a certain truncation that is within the o(h) margin. 
Although this correction is asymptotically small, it allows us to obtain several bounds in the 
convergence proof and also enables to construct transition probabilities that always remain 
in the unit interval (see ()3.12j) . below). So we now define 

(3.11) & k :=/3 k := 1 V fc ~ 1 ' ^ J , 1 < k < n, 



where 



and we set 




The reiterate once again that the function A n is chosen to ensure that the probabilities that 
are defined in (|3.12[) . below, remain in the unit interval. Although, this is clearly crucial for 
our analysis, in our numerical implementation we do not use this truncation and instead 
modify (|3.12[) to ensure that these are true probabilities. 

The above construction together with the conditional independence of the increments 
ensures the second moment matching. We now use the first order conditions (|3.10|) and 
p. 61) to construct the transition probabilities. Indeed, recall that by (|3.9I) . 

Pk :=P&(^ = 1) 

and rewrite (|3.10[) as 

p k exp U/hrj [(1 + &k) - otk-itk-i]) 

+(1 - p k ) exp (— s/hrj [(1 + a k ) + &k-i£k-i]j = exp(r/i). 
This implies that p k must be given by 

^ ^ _ exp(r/t + y/rjK afc_i£g_i) - exp(-y^ (1 + a k )) 

exp(y/rjh (1 + a k )) - cxp(- v / ?^i (1 + a k )) 
In view of the truncation introduced in (|3.11l) . p k G [0, 1] for all large n. 
We now recall that 

9k :=Pt ) 1 (4 y = l) 
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and use (13.61) to arrive at 



2 2(l + a k y k 2^(1 -p 2 ) (l + a k ) 
Since q k must take values in the unit interval, we modify it in the following way, 

1 , &k-l &Y V^M^-W-l) 1\ + 



(3.13) q k = I min i 1, - + ± . ^-1 



2 2(l + d fe ) sfc 2^(l-p 2 )(l + «fc) j / 
Set 

Then, we claim that S is a Markov process. Indeed, recall that the independence of the 
Brownian motions in (|2.1j) implies the conditional independence of the increments £ x and 
£ Y . Hence, 

(3.14) PW (tf = a,tf = 6 | S fc _i) = Pj&(# - «) P&ttf = 6). 

Moreover, in view of (13.11) and (13.2[) . the set 

{X<"> = 41 + c, F fc W = y W + d, = a, = 6} 



is empty unless c = 077/1 and d — br\h^f\ — p 2 , and in this case it is equal to {£ x — a, = 
b}. Therefore, the transition probabilities of the process S are determined by 

(# = l,# = l|S*_i)= Pfc g fc . 

Moreover, there is a simple transformation between and 

-ft •— , Y k ,t, k ,t, k )■ 

Hence, one may consider the process S as the basic approximating Markov process. 



4. Main convergence result 

In this section, we first briefly recall the concept of weak convergence of probability 
measures and then state our main convergence result. For more information on weak con- 
vergence, we refer the reader to the books of Billingsley [3] and Ethier & Kurtz [T4"] . 

For any cadlag stochastic process {Z(t)}J_ with values in some Euclidean space M. d , 
let F z be the distribution of Z on the canonical space O([0,T];R d ) equipped with the 
Skorohod topology (for details see [3]), i.e. for any Borel set D C D([0,T];R d ), V Z (D) = 
P{Z G D}. For a sequence of R d - valued, stochastic processes 

Z {n) 

we use the notation 

Z( n ) =>■ z to indicate that the probability measures P zl "' , converge vaguely to P z on the 
space B>([0,T];R d ). 

We are now ready to state the main convergence theorem which is the main theoretical 
foundation of our numerical scheme. It will be proved in Section [6] 

Theorem 4.1. For any n € N, let p( n ) be the probability measure defined by (|3.14l) . Con- 
sider the stochastic processes {4«t/T]}?=Q> { X [nt/T]}t=0 anc ^ i^[nt/T] }f=o under p( n ). Let 
(a;, y) be the unique solution of (|2.1l) . Then, 

(4-i) {(^/ti^S/t])}^ => {(*.,»*)}L> 

and 

(4-2) {(4«"/T]' r £/T])}Lo =* {(xt,yt)}to 

on the space D([0,T]) x D([Q,T]). 
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Remark 4.2. We emphasize that our approximation method using correlated random walks 
and the above convergence result can easily be extended to more general multidimensional 
diffusions. The key idea is the introduction of X type processes which differ from the original 
random walk X only by a predictable process a times the increment £ . We then use this 
freedom (namely, the function a) to construct transition probabilities that match the first 
and the second conditional moments of the original diffusion. The approximating process 
has the essentially same dimension as the original diffusion process. However, we need to 
augment the state space by adding the increments like £ x . But these increments take values 
in the discrete set { — 1, +1} so does not increase the complexity of the approximation. □ 

Our next remark is towards American options. 

Remark 4.3. In general, the usual weak convergence is not sufficient for the convergence 
of American options prices. Indeed, the latter also require the "good" behavior of the 
nitrations. In his unpublished manuscript (see pQ, Sections 15-16), David Aldous introduced 
the concepts of extended weak convergence to address this problem. Briefly his definition 
is as follows. A sequence : tt n -> D([0,T];R d ), extended weak converges to a stochastic 
process Z : Q — > D([0, T]; M d ), if for any k and continuous bounded functions ...,tpk & 
C(B([0,T];K rf )), 

(Z^ n \Z n ' 1 ,...,Z n ' k ) {Z,Z^,...,Z^) in D([0,T];K d+fe ), 

where for any t<T,l<i<k and n 6 N, 

= E^(MZ (n) ) I Jf '"'), zf = E(MZ) I 7f ), 

E^ n ' denotes the expectation on the probability space on which Z^ is defined and E 
denotes the expectation on the probability space on which Z is defined. In the formulas 
above F z{n) and T z are the filtration which are generated by Z^ and Z, respectively. The 
notion of extended weak convergence provides (in addition to the standard weak convergence 
of stochastic processes) convergence of nitrations. In particular, Aldous proved (see pQ 
Section 17) that under uniform integrability of the payoffs, extended weak convergence 
implies convergence of optimal stopping values. However, it is known that when the proof 
of weak convergence relies on martingale techniques (like our proof), then the standard 
weak convergence implies extended weak convergence. For details, we refer the reader to 
PQ, Section 21. □ 

5. Discrete pricing equations 

In this section, we apply the approximation developed in Section [3] to price American 
put and lookback options. 

5.1. American Put. Consider an American Put option with a strike price K. We are 
interested in approximating its value given by, 

V= sup E(e- rT (K - S T ) + ) , 

i"G7]o,T] 

where 7[o,t] is the set of all stopping times with respect to the filtration generated by 5, 
with values in the set [0, T]. We approximate the discounted stock price by the the discrete 
time martingales 

- ( n \ n 

{e- rkh e x * } fc=0 , nel, 
constructed in Section[3] For any n £ N, let T n be the set of all stopping times with respect 
to the filtration Fk (again constructed in Section [3]), with values in the set {0, 1, n}. 
Define, 

V (n) := maxE("' f e - rTh {K - S e x ^)+) . 
In view of Theorem 14.11 and Remark 14.31 we directly conclude that 

lim V {n) = V. 
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Next, we describe a dynamical programming algorithm for the calculation of V^ n \ Ob- 



serve that for a given k £ {0, . . . , n} the random variables and take values on the 



grid 

x + (21 - k) y/rjh, < Z < fc, 
y + (2m - k)y/r)(l - p 2 )h, < m < k, 
respectively. For non-negative integers m,l < k < n and £, Xl £,y € {~ 1, let 

be the value of the option at time k when the Markov process is given by 

E k = (xj?\Y£ n \&,&)=F h (!,m,Z x ,ty) 

:= (x Q + (21 - k)y^jh, yo + (2m - k)^ V (l - p 2 )h, £„). 

The above function i 7 ^. is invertible with an inverse FjT . We sometimes, with an abuse of 
notation, write 

for any four tuple S given by Ffc_i (/, to, ^) for some (i, to, 6, £ y ). With this convention, 
it is not straightforward to state the dynamic programming equation (see, for instance |26] , 
Chapter 1), 



— fe-i — — 



(5.1) v£>(5) = max | (k - S exp(l fc _ 1 )) + , E (n) [v£ n) {2 k ) 

We continue by rewriting the dynamic programming equation in an algorithmic manner. 
In view of (|3.11l) - (|3.13j) . for any 1 < k < n and < I, to < k — 1, we define, 



X k := x + (2l-kWr]h, 

y k := y + (2m-k)y/r)(l- p 2 )h, 

where both of the above are functions of (l,m), but this dependence is supressed in the 
notation. Similarly, we define two probabilities 

exp(rh + <JrjK # fc _i 60 - exp(-y/rjh * fc ) 



Pk(l,m,£ x ,£ y ) := 



exp(y/r]h * fc ) - exp(-V^ *k) 

1 Q!fe_l(i - £ x , TO - £}, Vh Py,k 



q k (l,m,£ x ,Q := [minil, , _^ 



where = and 



a k (l,m) 



max(A„,cr 2 (X k -i,y k -i)) - 1 



2 



:= l + a k n, (l,m), 

As remarked earlier, in our actual numerical codes, we simply define a — (a 2 — l)/2 without 
the truncation with A n and instead truncate p k , above, to ensure that it says within the 
unit interval. 
Observe that 

Pk (l,m,U^y) = lP (n) (Cf = 1 I Sfc_i = F k -i(l,m,i,c,tv) ) , 
q k (l,m,Z x ,Z y ) = P<">(# = l|S fc _;t = F fc _i(i,m, 

Moreover, 

Pfc-i (tf = 1 . # = 1 ) =p fc a,m,6,^) ffk(i,m, 
One can easily obtain expressions for other three probabilities as well. 
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We now are ready to restate the dynamic programming equation (I5.ip . Indeed, (I, m, £ x , £ y ) 
is the unique solution of the following recursive relations, 

V,i n) [I, m, £ X) S s ) = (k-S exp (yfijh (X n + ) + , 

and for 1 < k < n. 



(I, m, £ x , £ y ) = max | (k - 5 exp (^Jrjh {X k - 



where 

£(T/ fe (n) ) - EW [vf } (S fc ) I S fc _! = F k ^{l t m,^ v ) 
1 

= E P £-i = 2* - 1, 4 F = 2j - 1) V k {n) (l + i,m + j, 2i — 1, 2j - 1) 

i,j=0 

l 

= ^ [1 - t + (2t- l)j>it(i,ro, £»,£„)] [1 - j + (2j - l)g fc (/, m, 
»>j=o 

x V k [n \l + i,m + j,2i- l,2j - 1). 
Then, our approximation is simply given by, 

K = ^ (n) (0,0,0,0). 

5.2. Lookback options. Consider a lookback put option with a fixed strike K, i.e., an 
option with payoff (K — mino<«T <5t) + ■ Again, we want to approximate the price 

V = E | e" rT fiT- min ^ 

\ V 0<t<T 

Since the running minimum of the processes 

{exp^)}^, neN, 

lies on a grid, we will use these processes instead of the martingale exp(A > ^ n ' ) . The advantage 

of the processes exp(-X^ ) becomes clear, when we describe the dynamical programming 
algorithm, below. 
We set, 

(5.2) V-W = (e- rT (k - S exp (^^jj . 

By Theorem 14. II we conclude that yw converges to V'. 
First, we observe that the random variable 

i 

z k := min V^f 

~ - 3=1 

takes values on the grid {—k, 1 — k, 0}. 

Using the notations and the conventions of the previous subsection, for < k < n, we 
let V/. (I, m, z, £, x ,£y) be the option price at time k. The extra state variable z denotes the 
value of the running minimum z k at time k. Then, V"' n ) is the unique solution of 

Vt ] (l,m,z,£ x ,Z y ) = (K-S exp(-^z)y 

and for 1 < k < n, 

V k ( ™\(l, TO > z > ^' £v) = max I ( A ' _ exp (-y/rjh z)) , £ (F fc (n) ) 1 , 
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where 

£ (vt } ) = E P& (a x = a - 1, = sy - 1) 

i,j=0 

x F fc (n) (/ + i,m + j, z + X{i=o,z+2i=k-i},2i - l,2j - l) , 
and xq is the characteristic set of Q. Finally, 

v n = v {n) (0,0, 0,0,0). 
6. Proof of Theorem 2.1 



In this section we provide a proof of Theorem 14.11 Our main tool is the martingale 
convergence result of Theorem 7.4.1 in [T4] . 

In view of (|3.1[) - (|3.4[) and (I3.11[) . we have the following inequality for all sufficiently large 

n, 

i^ n) i>i4 B) i-5(i4 B) i + in (B) i + i). 



Therefore, 



in (n) i>in (n) i-i(i4 n) i + i^ (n) i+D- 



(6.1) \X^\ + \Y^\<3(\X^\ + \Y^\ + 1), fc = 0,l,...,n. 

This together with (|3.3|l - (|3.4|l and p.llj) implies that there exists a constant c > 
satisfying, 

(6.2) \X™ - Xi n) \ + \Yl n) - y fc (n) ] < 0(1 + |Yfc( " )|} , fc = 0,l,...,n. 



We conclude that in order to prove Theorem 14.11 it sufficient to establish that 

(6-3) {(x}:l T] ,Y{ n 1 /T] )}Lo =* 

For any < fe < n, set 

j=i j'=i 

M „,x = £<«) _ ^ M „, y = y(n) _ A n,y^ 

k k 

A n,*, x = ^2 E (»\ ((m;-* - m;^) 2 ) , A^ ,y,!/ = £ 4-i (( M "' y - M "-i) 2 ) . 

3=1 3=1 
fe 

^T' 9 = E E !"-i (( M ?" - M "- x ( m j' v - m j~i)) ■ 

i=i 

Notice that the processes A n,x , A n ' y , A n - X ' x . A n ' V:V , A n ' x,v are predictable and the processes 
M n,x ,M n ' y are martingales. 

We now fix a large N > and define the stopping times by, 

a n = min{fe : \X ( k n) \ + |f fe (n) | > iV} A n, n e N. 

Using (|3.1[) . fl3.2[) and (|6.2[) , we conclude that for all k < a n , 

X ( k ] ~ = 0(V V^), and ff> - f£> = 0(l/y/n), 

where in this section o(-) and O(-) are defined uniformly in space, i.e., 0(1/ \/n) is a function 
which is bounded by a constant times y 7 ?! and i/n o(l/i/n) converges uniformly to zero as 
n tends to infinity. 
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By Theorem 7.4.1 in [14], (|6.3j) would result from the following relations, 

fc-i 



(6.4) 

(6.5) 

(6.6) 

(6.7) 
(6.8) 



lim max 

n^-oo l<k<<T n 



lim max 

n— yoo l<k<<7 n 



lim max 

n^-oo l<k<<T n 



lim max 

n^-oo l<k<<T n 



i=0 



i=0 



fc-1 



a.s., 
a.s., 
= a.s., 



i=0 



fe-1 



^(1-p 2 )/* X) ff2 W n) ^ (n) ) 



i=0 



a.s. 



lim max | J 4^' :r ' !/ | = a.s. 



n— too l<k<a n 

The rest of the proof is devoted to the verification of the above identities. 
We start with a proof of (|6.4p . Since a 2 (x,y) is Lipschitz continuous, (|3.1j) . (|3.2p and 
p. lip implies that 

|a*-a*-i)|=0(Vfc). 
In view of (|6.1|) . for k < a n , we have, 

~2 <"fe <c(JV+l) 
for some constant c. Since the event k < a n is J-fc_i-measurable, 

PW (4 X = 1 and k < a n | S fc _i) = X {fe<CTn} P (n) fe X - 1 I 3fe-i) = X { /c< CT „} Pk- 

We now use the above estimates, the definition (|3.12[> of the transition probability p k and 
Taylor expansion. Then, on the set k < a n , 

= 1 I 5fe-i) 

_ rfe + (1 + a k -i^k-i + ~ V h (V 2 + Q^fc) + °( h ) 

(1 + a k ) + o(h) 
_ rh + V^fe (1 + a fc -i^_i + &fc) - ^ (1/2 + a k ) 
2^jh(l + a k ) 

(6-9) = - + , . . 4-i + on ~ . \ + o(/i). 

2 2(1 + a*,) 2(l + afe) 

We thus conclude that on the event k < <r n , the following estimate holds, 

= E k n \ [(1 + a k )$ - & k -xg-.i] 



e 3 i 



(1 + a k ) (2P(")(^ = 1 I S fc _ x ) - 1) - i-i4 X i 
= rh-T](- +&k)h + o(h) 
= f , x (xi%Y£\)h + o(h), 

where the last equality follows from the definition of a, the Lipschitz continuity of (J,(x, y)and 
(16.21) . Then, (|6.4p follows directly from the above estimate. 

We continue with a proof of (I6.5p . We start with the definition of q k and use the 
truncation introduced in (|3.1ip . On k < <r n , this fields the following estimate, 



2PW (4 Y = 1 I 3fe-i) - 1 = ff-kJ : 

1 + Otk 
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As before we directly estimate the on k — 1 < a n , 

= ^v(i-P 2 )h ((1 + d fc ) (2PW (4 y = 1 1 s fc _ x ) - l) - d fc Ci) 



Again, the last equality follows from (|6.2[) and the fact that p y (x, y) is Lipschitz continuous. 
This completes the proof of (|6.5p . 

We continue with the quadratic estimates. Indeed, by (|6.9|) . on k < a n , 

2PW (tf = llSfc-x) - 1 = f^— ^ , + o(^). 
Since A n,a; is predictable, on k < a n , 

4" } i ((JiC" - <-i) 2 ) = 4- } i ((4 n) - *&) a ) + o(h) 

and 

= V h ((1 + a k ) 2 + (a k ^) 2 - 2d fc _i(l + (2P (n) (£? = l|H fc _i) - 1)) 

= 77/1 (i + 2d4 n) ) 

and (|6.6|) follows. The relation (|6.7|) is is proved similarly. 

It remains to establish (|6.8|) . The processes A™'*, J 4™' 1 ' are predictable. Thus, from (|3.14[) 
it follows that, on k < a n , 

= E& ((if 1 - *&) (yW _ yW)) + o(/l) 

= E W - XW )E& (if > - Yt\) + o(h) 

= o(h). 

□ 

7. Numerical Results 

In this section, we present numerical results from our model for European and American 
Vanilla, Lookback, Geometric and Arithmetic Asian options under the Heston dynamics. 
Our computations are obtained by a direct implementation of the methodology described 
in the previous sections. In particular, we explicitly refrained from using known numerical 
techniques that improve the performance of the trees. This is done to ensure the replicability 
of our reported results. We will also provide our codes upon request. 

7.1. Vanillas. In Tables 1, 2, and 3 we use the same parameter sets as in Beliaeva and 
Nawalkha [2], i.e. for European call and put options: strike K = 100; initial stock prices: 
So = 90, 95, 100, 105, 110, maturities: T = 1 month, 3 months, and 6 months; initial volatil- 
ity values: ^Jva = 0.2,0.3,0.4; interest rate: r — 0.05; vol of vol: i] — 0.1; mean reversion 
rate: k — 3; long run vol: 9 = 0.04; and correlation: p — —0.7. For American put op- 
tions: K = 100, S = 90, 100, 110; T = 1 month, 3 months, and 6 months; ^ = 0.2,0.4; 
p = -0.1, -0.7; r = 0.05; 77 = 0.1; K = 3, 9 = 0.04. 

Tables 1 and 2 show the convergence of European put and call prices computed by our 
method compared to the closed form solutions of Heston [T7] . In the European case, one can 
calculate errors as Heston's solution is available in closed form. The option prices computed 
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for the number of time steps N = 200, 350 and 500 illustrate very good convergence to the 
closed form solution as reported in Tables 1 and 2. Furthermore, one verify that the put-call 
parity holds exactly for option prices at each of these time steps sizes. Clearly, this is the 
outcome of the fact that our price process in any step size is a martingale. 

Table 3 reports the difference between the American put prices obtained from our method 
and those obtained by the Control Variate (CV) technique of [2J. The table shows that our 
numbers are in good agreement with those obtained by the CV method. The first three 
largest differences between the models are (0.27%, 0.26%, 0.22%) and on average there is a 
difference of 0.10% per option. We should point out to the reader that the CV technique 
computes the value of the put option modifies the computed prices via the fromula 
CV American Price = Tree American + (Closed Form Euro - Tree Euro). 

According to [2], this method is particularly useful for longer maturity options. 

Chockalingam and Muthuraman [B] develop a partial differential equations (PDE) based 
finite difference method to price American options under stochastic volatility. More specif- 
ically, they transform the free boundary problem resulting from the pricing of American 
options into a sequence of fixed-boundary problems of European type. The prices listed in 
Table 4 and 5 are taken from [6] as a benchmark for our tree based method. The authors 
provide the values arising from the projected successive over relaxation (PSOR) method 
and the component-wise splitting (CS) method. They state that other PDE based methods 
(see Ikonen and Toivanen [2U] for a detailed analysis) fall between these two in terms of 
speed/accuracy and ease of implementation. As test parameters, they use the most common 
parameter values for American options under the Heston dynamics in the PDE-based liter- 
ature: K = 10, r = 0.1, T) = 0.9, k = 5.0, 9 = 0.16, and p = 0.1, T — 0.25, ^ = 0.25, 0.5. 
Following 6], we take the prices computed by Ikonen and Toivanen 20 (using the CS 
method together with a very fine grid) as the reference values. From Table 4 and 5, one can 
clearly conclude that our results for both N — 250 and N = 350 are very close to reference 
values. 

7.2. Exotics. Our numerical experimentation confirms that backward recursion yields quite 
fast and accurate results for the two dimensional problems like European and American 
Vanilla option pricing problems. However, our numerical experimentation also reveals that 
the straightforward application of the recursive method takes too long on a personal com- 
puter when another continuous variable is introduced to price an exotic option. Hence, 
in order to substantially speed up the computations, we use our discrete equations as a 
discretization scheme for our Monte Carlo (MC) simulation. In other words, we carry out 
the MC simulation on the tree. 

It is also important to note that our main concern in this section is to show the pure 
application of our computation method. There are many well known techniques in the 
literature which improve the speed and the accuracy of tree and MC methods. However, as 
in the backward recursion we refrain using any of these techniques. 

Below we outline results for the geometric, arithmetic Asian and for lookback options. 

We start with the geometric Asian and let 



be the geometric mean of St over time t during [0,T]. Then the payoff of a fixed strike 
geometric Asian call is given by max(GT — K,0). Kim and Wee [21] provide semi-closed 
solutions for the price of geometric Asian options under the Heston model. We compare 
our results with theirs. 

Table 5 displays a comparison between prices from the semi-closed solution and those 
from our MC simulation on tree with TV = 300 and number of simulations (NumSim)= 
10 5 , 5 * 10 5 , 10 6 . As benchmark prices, we use the values given in Table 5 from [2T] for the 
parameter values: S = 100, u Q = 0.09, r = 0.05, k = 1.15, 6 = 0.348, p = -0.64, 77 = 0.39. 



G T = exp - 
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As it is clear from the table, our numerical scheme provides a very good approximation 
for the analytical prices. For NumSim = 10 6 , we get three largest percentage errors 
as (0.40%, 0.34%, 0.28%) and average percentage error is 0.11%. Table 8 shows the 95% 
confidence intervals for the prices computed for different numbers of simulations. 

Table 8 includes our results for arithmetic Asian options under the Heston model. We 
carry out the simulations as in the same way described previously. Let 



be the arithmetic average of St over time t during [0,T]. Then the payoff of a fixed strike 
arithmetic Asian call is given by max(AT — K, 0). Pages and Printems |25] use the functional 
quantization based quadrature formula to price vanilla Calls and Asian Calls in the Heston 
model. The numbers computed from MC method, Romberg log-extrapolation, and K- 
interpolation of Romberg and their standard deviations in the parenthesis are tabulated for 
comparison (see Table 4 in [25^ for a more detailed explanation of the results). We test 
our model using the numbers reported in their paper. As one can observe from Table 8, 
our prices together with the confidence intervals are in accordance with the only reference 
values for arithmetic Asian options under the Heston dynamics which can be found in the 
literature. 

It is clear that when we price a lookback option using backward recursion, we also need 
another continuous variable holding the running max or min. But in this case, we can 
constrain this variable to take values on a tree as well. However, it still remains more 
efficient to apply our MC method on the tree. Table 9 presents numerical results obtained 
by the standard MC method and our numerical method for fixed strike lookback Call 
options. As comparison we used simple Monte Carlo simulations based on a Euler method. 
The table contains prices for N — 3000 and NumSim — 10 5 . As one can see from the 
last column, the numbers obtained from our numerical method differ only slightly from the 
prices computed by the Euler MC method. 

Table 10 provides average running times for the options in Tables 1-9. The computer 
used is a standard laptop with an Intel Core i7 M620@2.67 GHz CPU and a 4-GB memory. 
The algorithm was implemented in MATLAB. 



In this paper, we have developed a recombining tree approximation of the Heston model. 
Our approach is very general and applies to all stochastic volatility models with a factor 
equations. For low dimensional European and the American option equation can be solved 
by a straightforward backward recursion. 

We have done extensive numerical experimentation with the resulting pricing equations. 
These results, reported in the previous section, confirms the efficiency of the method. 
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Table 1 . Convergence of European put prices versus analytical solution 
of Heston [17] . Parameters: K = 100, r = 0.05, r\ = 0.1, K = 3.0, 9 = 0.04, 
and p = -0.7. 



S(0) 


V "o 


T 


Tree 


Analytical 
Solution 


Error % 


N = 200 


N = 350 


N = 500 


N = 200 


N = 350 


N = 500 


90 


0.2 


0.0833 


9.6541 


9.6533 


9.6533 


9.6533 


0.01 


0.00 


0.00 


95 


0.2 


0.0833 


5.2059 


5.2084 


5.2077 


5.2074 


-0.03 


0.02 


0.01 


100 


0.2 


0.0833 


2.0953 


2.0960 


2.0965 


2.0971 


-0.08 


-0.05 


-0.03 


105 


0.2 


0.0833 


0.6082 


0.6047 


0.6050 


0.6053 


0.48 


-0.10 


-0.06 


110 


0.2 


0.0833 


0.1267 


0.1271 


0.1270 


0.1265 


0.11 


0.48 


0.35 


90 


0.3 


0.0833 


9.9913 


9.9900 


9.9900 


9.9905 


0.01 


0.00 


0.00 


95 


0.3 


0.0833 


6.0147 


6.0170 


6.0162 


6.0155 


-0.01 


0.02 


0.01 


100 


0.3 


0.0833 


3.1308 


3.1288 


3.1290 


3.1302 


0.02 


-0.05 


-0.04 


105 


0.3 


0.0833 


1.4001 


1.3955 


1.3955 


1.3967 


0.25 


-0.08 


-0.09 


110 


0.3 


0.0833 


0.5365 


0.5374 


0.5372 


0.5367 


-0.05 


0.13 


0.09 


90 


0.4 


0.0833 


10.5687 


10.5670 


10.5668 


10.5668 


0.02 


0.00 


0.00 


95 


0.4 


0.0833 


6.9357 


6.9363 


6.9352 


6.9335 


0.03 


0.04 


0.02 


100 


0.4 


0.0833 


4.1893 


4.1864 


4.1861 


4.1852 


0.10 


0.03 


0.02 


105 


0.4 


0.0833 


2.3280 


2.3232 


2.3229 


2.3222 


0.25 


0.04 


0.03 


110 


0.4 


0.0833 


1.1893 


1.1897 


1.1893 


1.1882 


0.09 


0.13 


0.09 


90 


0.2 


0.25 


9.5736 


9.5693 


9.5694 


9.5698 


0.04 


0.00 


0.00 


95 


0.2 


0.25 


5.9691 


5.9685 


5.9693 


5.9692 


0.00 


-0.01 


0.00 


100 


0.2 


0.25 


3.3742 


3.3774 


3.3794 


3.3770 


-0.08 


0.01 


0.07 


105 


0.2 


0.25 


1.7420 


1.7393 


1.7402 


1.7410 


0.06 


-0.10 


-0.05 


110 


0.2 


0.25 


0.8290 


0.8249 


0.8253 


0.8259 


0.37 


-0.13 


-0.08 


90 


0.3 


0.25 


10.5941 


10.5879 


10.5882 


10.5893 


0.04 


-0.01 


-0.01 


95 


0.3 


0.25 


7.3343 


7.3327 


7.3329 


7.3316 


0.04 


0.02 


0.02 


100 


0.3 


0.25 


4.8279 


4.8331 


4.8340 


4.8310 


-0.06 


0.04 


0.06 


105 


0.3 


0.25 


3.0420 


3.0379 


3.0391 


3.0388 


0.11 


-0.03 


0.01 


110 


0.3 


0.25 


1.8368 


1.8320 


1.8319 


1.8325 


0.23 


-0.03 


-0.03 


90 


0.4 


0.25 


11.8375 


11.8281 


11.8288 


11.8287 


0.07 


0.00 


0.00 


95 


0.4 


0.25 


8.8120 


8.8081 


8.8070 


8.8035 


0.10 


0.05 


0.04 


100 


0.4 


0.25 


6.3762 


6.3790 


6.3786 


6.3735 


0.04 


0.09 


0.08 


105 


0.4 


0.25 


4.5066 


4.5005 


4.5004 


4.4976 


0.20 


0.06 


0.06 


110 


0.4 


0.25 


3.1099 


3.1035 


3.1025 


3.1011 


0.28 


0.08 


0.05 


90 


0.2 


0.5 


9.7547 


9.7545 


9.7606 


9.7572 


-0.03 


-0.03 


0.04 


95 


0.2 


0.5 


6.7258 


6.7248 


6.7185 


6.7199 


0.09 


0.07 


-0.02 


100 


0.2 


0.5 


4.4355 


4.4369 


4.4320 


4.4312 


0.10 


0.13 


0.02 


105 


0.2 


0.5 


2.8077 


2.8159 


2.8100 


2.8107 


-0.11 


0.18 


-0.02 


110 


0.2 


0.5 


1.7286 


1.7289 


1.7275 


1.7240 


0.27 


0.28 


0.20 


90 


0.3 


0.5 


11.0786 


11.0792 


11.0845 


11.0807 


-0.02 


-0.01 


0.03 


95 


0.3 


0.5 


8.2445 


8.2422 


8.2367 


8.2363 


0.10 


0.07 


0.00 


100 


0.3 


0.5 


5.9835 


5.9830 


5.9784 


5.9763 


0.12 


0.11 


0.04 


105 


0.3 


0.5 


4.2450 


4.2504 


4.2449 


4.2443 


0.02 


0.15 


0.02 


110 


0.3 


0.5 


2.9647 


2.9640 


2.9623 


2.9582 


0.22 


0.20 


0.14 


90 


0.4 


0.5 


12.6195 


12.6199 


12.6231 


12.6171 


0.02 


0.02 


0.05 


95 


0.4 


0.5 


9.9373 


9.9318 


9.9260 


9.9223 


0.15 


0.10 


0.04 


100 


0.4 


0.5 


7.7110 


7.7069 


7.7017 


7.6965 


0.19 


0.13 


0.07 


105 


0.4 


0.5 


5.9065 


5.9075 


5.9015 


5.8978 


0.15 


0.17 


0.06 


110 


0.4 


0.5 


4.4841 


4.4806 


4.4779 


4.4716 


0.28 


0.20 


0.14 
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Table 2. Convergence of European call prices versus analytical solution 
of Heston Q7|. Parameters: K = 100, r = 0.05, 77 = 0.1, K = 3.0, 6 = 0.04, 
and p — —0.7. 



S(0) 


V "0 


T 


Tree 


Analytical 
Solution 


Error % 


N = 200 


N = 350 


N = 500 


N = 200 


N = 350 


N = 500 


90 


0.2 


0.0833 


0.0699 


0.0691 


0.0691 


0.0691 


1.13 


-0.05 


0.02 


95 


0.2 


0.0833 


0.6217 


0.6242 


0.6235 


0.6232 


-0.23 


0.17 


0.06 


100 


0.2 


0.0833 


2.5111 


2.5118 


2.5122 


2.5129 


-0.07 


-0.04 


-0.02 


105 


0.2 


0.0833 


6.0240 


6.0205 


6.0208 


6.0211 


0.05 


-0.01 


-0.01 


110 


0.2 


0.0833 


10.5425 


10.5429 


10.5428 


10.5423 


0.00 


0.01 


0.00 


90 


0.3 


0.0833 


0.4071 


0.4058 


0.4058 


0.4063 


0.20 


-0.12 


-0.12 


95 


0.3 


0.0833 


1.4305 


1.4328 


1.4320 


1.4313 


-0.06 


0.10 


0.05 


100 


0.3 


0.0833 


3.5466 


3.5446 


3.5448 


3.5460 


0.02 


-0.04 


-0.04 


105 


0.3 


0.0833 


6.8159 


6.8113 


6.8113 


6.8125 


0.05 


-0.02 


-0.02 


110 


0.3 


0.0833 


10.9523 


10.9532 


10.9530 


10.9525 


0.00 


0.01 


0.00 


90 


0.4 


0.0833 


0.9845 


0.9828 


0.9826 


0.9826 


0.19 


0.02 


0.00 


95 


0.4 


0.0833 


2.3515 


2.3521 


2.3510 


2.3493 


0.10 


0.12 


0.07 


100 


0.4 


0.0833 


4.6051 


4.6022 


4.6019 


4.6010 


0.09 


0.03 


0.02 


105 


0.4 


0.0833 


7.7438 


7.7390 


7.7387 


7.7380 


0.08 


0.01 


0.01 


110 


0.4 


0.0833 


11.6051 


11.6055 


11.6051 


11.6040 


0.01 


0.01 


0.01 


90 


0.2 


0.25 


0.8158 


0.8115 


0.8116 


0.8120 


0.47 


-0.05 


-0.05 


95 


0.2 


0.25 


2.2113 


2.2107 


2.2116 


2.2114 


-0.01 


-0.03 


0.01 


100 


0.2 


0.25 


4.6164 


4.6196 


4.6216 


4.6192 


-0.06 


0.01 


0.05 


105 


0.2 


0.25 


7.9842 


7.9815 


7.9824 


7.9832 


0.01 


-0.02 


-0.01 


110 


0.2 


0.25 


12.0712 


12.0671 


12.0675 


12.0682 


0.03 


-0.01 


-0.01 


90 


0.3 


0.25 


1.8363 


1.8301 


1.8305 


1.8316 


0.26 


-0.08 


-0.06 


95 


0.3 


0.25 


3.5766 


3.5750 


3.5751 


3.5738 


0.08 


0.03 


0.03 


100 


0.3 


0.25 


6.0701 


6.0753 


6.0762 


6.0732 


-0.05 


0.03 


0.05 


105 


0.3 


0.25 


9.2842 


9.2802 


9.2813 


9.2810 


0.04 


-0.01 


0.00 


110 


0.3 


0.25 


13.0790 


13.0742 


13.0741 


13.0747 


0.03 


0.00 


0.00 


90 


0.4 


0.25 


3.0797 


3.0703 


3.0710 


3.0709 


0.29 


-0.02 


0.00 


95 


0.4 


0.25 


5.0542 


5.0503 


5.0493 


5.0457 


0.17 


0.09 


0.07 


100 


0.4 


0.25 


7.6184 


7.6212 


7.6208 


7.6157 


0.04 


0.07 


0.07 


105 


0.4 


0.25 


10.7488 


10.7428 


10.7426 


10.7399 


0.08 


0.03 


0.03 


110 


0.4 


0.25 


14.3521 


14.3457 


14.3447 


14.3433 


0.06 


0.02 


0.01 


90 


0.2 


0.5 


2.2237 


2.2235 


2.2296 


2.2262 


-0.11 


-0.12 


0.15 


95 


0.2 


0.5 


4.1948 


4.1938 


4.1875 


4.1889 


0.14 


0.12 


-0.03 


100 


0.2 


0.5 


6.9045 


6.9060 


6.9010 


6.9002 


0.06 


0.08 


0.01 


105 


0.2 


0.5 


10.2767 


10.2849 


10.2790 


10.2797 


-0.03 


0.05 


-0.01 


110 


0.2 


0.5 


14.1976 


14.1979 


14.1965 


14.1930 


0.03 


0.03 


0.02 


90 


0.3 


0.5 


3.5476 


3.5483 


3.5535 


3.5497 


-0.06 


-0.04 


0.11 


95 


0.3 


0.5 


5.7135 


5.7112 


5.7057 


5.7053 


0.14 


0.10 


0.01 


100 


0.3 


0.5 


8.4525 


8.4520 


8.4474 


8.4453 


0.09 


0.08 


0.03 


105 


0.3 


0.5 


11.7140 


11.7194 


11.7140 


11.7133 


0.01 


0.05 


0.01 


110 


0.3 


0.5 


15.4337 


15.4330 


15.4313 


15.4272 


0.04 


0.04 


0.03 


90 


0.4 


0.5 


5.0885 


5.0889 


5.0921 


5.0861 


0.05 


0.06 


0.12 


95 


0.4 


0.5 


7.4063 


7.4008 


7.3950 


7.3913 


0.20 


0.13 


0.05 


100 


0.4 


0.5 


10.1800 


10.1759 


10.1707 


10.1655 


0.14 


0.10 


0.05 


105 


0.4 


0.5 


13.3755 


13.3765 


13.3705 


13.3668 


0.07 


0.07 


0.03 


110 


0.4 


0.5 


16.9532 


16.9496 


16.9469 


16.9406 


0.07 


0.05 


0.04 
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Table 3. Comparison of American put prices calculated with our method 
and with the control variate technique of Beliaeva and Nawalkha [5] . Pa- 
rameters: K = 100, r = 0.05, tj = 0.1, K = 3.0, = 0.04, and p = -0.7. 



5(0) 


P 




T 


Tree 


Control Variate 


Difference % 


N = 250 


N = 200 


90 


-0.1 


0.2 


0.0833 


10.0000 


10.0000 


0.00 


100 


-0.1 


0.2 


0.0833 


2.1236 


2.1254 


-0.08 


110 


-0.1 


0.2 


0.0833 


0.1090 


0.1091 


-0.05 


90 


-0.7 


0.2 


0.0833 


10.0000 


9.9997 


0.00 


100 


-0.7 


0.2 


0.0833 


2.1249 


2.1267 


-0.08 


110 


-0.7 


0.2 


0.0833 


0.1273 


0.1274 


-0.07 


90 


-0.1 


0.4 


0.0833 


10.7123 


10.7100 


0.02 


100 


-0.1 


0.4 


0.0833 


4.2194 


4.2158 


0.08 


110 


-0.1 


0.4 


0.0833 


1.1666 


1.1667 


-0.01 


90 


-0.7 


0.4 


0.0833 


10.6843 


10.6804 


0.04 


100 


-0.7 


0.4 


0.0833 


4.2183 


4.2140 


0.10 


110 


-0.7 


0.4 


0.0833 


1.1942 


1.1939 


0.02 


90 


-0.1 


0.2 


0.25 


10.1713 


10.1706 


0.01 


100 


-0.1 


0.2 


0.25 


3.4729 


3.4747 


-0.05 


110 


-0.1 


0.2 


0.25 


0.7726 


0.7736 


-0.13 


90 


-0.7 


0.2 


0.25 


10.1222 


10.1206 


0.02 


100 


-0.7 


0.2 


0.25 


3.4790 


3.4807 


-0.05 


110 


-0.7 


0.2 


0.25 


0.8405 


0.8416 


-0.13 


90 


-0.1 


0.4 


0.25 


12.1880 


12.1819 


0.05 


100 


-0.1 


0.4 


0.25 


6.5023 


6.4964 


0.09 


110 


-0.1 


0.4 


0.25 


3.0952 


3.0914 


0.12 


90 


-0.7 


0.4 


0.25 


12.1245 


12.1122 


0.10 


100 


-0.7 


0.4 


0.25 


6.4989 


6.4899 


0.14 


110 


-0.7 


0.4 


0.25 


3.1512 


3.1456 


0.18 


90 


-0.1 


0.2 


0.5 


10.6521 


10.6478 


0.04 


100 


-0.1 


0.2 


0.5 


4.6531 


4.6473 


0.12 


110 


-0.1 


0.2 


0.5 


1.6857 


1.6832 


0.15 


90 


-0.7 


0.2 


0.5 


10.5682 


10.5637 


0.04 


100 


-0.7 


0.2 


0.5 


4.6691 


4.6636 


0.12 


110 


-0.7 


0.2 


0.5 


1.7899 


1.7874 


0.14 


90 


-0.1 


0.4 


0.5 


13.3279 


13.3142 


0.10 


100 


-0.1 


0.4 


0.5 


8.0231 


8.0083 


0.18 


110 


-0.1 


0.4 


0.5 


4.5554 


4.5454 


0.22 


90 


-0.7 


0.4 


0.5 


13.2431 


13.2172 


0.20 


100 


-0.7 


0.4 


0.5 


8.0204 


7.9998 


0.26 


110 


-0.7 


0.4 


0.5 


4.6328 


4.6201 


0.27 
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Table 4. American put prices determined with our tree approach and 
finite difference methods. Parameters: K = 10, r = 0.1, rj = 0.9, k = 5.0, 
9 = 0.16, and p = 0.1, T = 0.25, y/vS = 0.25 







So 




VjrllU Olzie 


8 


9 


10 


11 


12 


PSOR 


(40,16,8) 
(60,32,66) 
(120,64,130) 
(240,128,258) 


2.0000 
2.0000 
2.0000 
o nnnn 

z.UUUU 


1.0952 
1.1037 
1.1064 

1 1 071 


0.4966 
0.5142 
0.5182 
U.oiyo 


0.2042 
0.2105 
0.2126 

U.ZLoo 


0.0838 
0.0815 
0.0819 
n neon 


splitting 


(60,32,66) 
(120,64,130) 
(240,128,258) 


2.0004 
2.0000 
2.0000 
9 noon 


1.1003 
1.1043 
1.1066 

±.±U 1 o 


0.4991 
0.5147 
0.5183 
u.o ly^ 


0.2035 
0.2104 
0.2126 
u.zioo 


0.0828 
0.0813 
0.0819 
u.uozu 


Trflri Qfnrm a tinn 

procedure 


f40 1 6 81 
(60,32,66) 
(120,64,130) 
(240,128,258) 


2.0000 
2.0000 
2.0000 
2.0000 


1.0952 
1.1035 
1.1063 
1.1071 


0.4966 
0.5142 
0.5181 
0.5193 


0.2042 
0.2105 
0.2126 
0.2133 


0.0838 
0.0815 
0.0819 
0.0820 


Our Tree Method 


N 














150 
250 
350 


2.0000 
2.0000 
2.0000 


1.1086 
1.1079 
1.1074 


0.5155 
0.5190 
0.5193 


0.2140 
0.2140 
0.2134 


0.0825 
0.0822 
0.0828 


Reference Value 




2.0000 


1.1076 


0.5200 


0.2137 


0.0820 
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Table 5. American put prices determined with our tree approach and 
finite difference methods. Parameters: K = 10, r = 0.1, r\ = 0.9, k = 5.0, 
= 0.16, and p = 0.1, T = 0.25, ^ = 0.5 







So 


IVlcLllUU. 


vjiicl oizie 


8 


9 


10 


11 


12 


PSOR 


(40,16,8) 
(60,32,66) 
(120,64,130) 
(240,128,258) 


2.0691 
2.0760 
2.0775 

0770 
LAM /y 


1.3139 
1.3292 
1.3320 

1 QQOO 

l.oozy 


0.7720 
0.7908 
0.7940 
o 70^1 


0.4293 
0.4442 
0.4467 

O A A IfK 
U.44 / 


0.2324 
0.2405 
0.2419 

O O/IO/I 
U.z4z4 


I Am Y~\ An ATI Til 71 (J A 

splitting 


(AO 1 A SJ^i 

(60,32,66) 
(120,64,130) 
(240,128,258) 


2.0676 
2.0758 
2.0774 

Z.U ( OU 


1.3094 
1.3287 
1.3317 

1 


0.7646 
0.7900 
0.7936 

O 7Q/1Q 

u. ( y^y 


0.4232 
0.4435 
0.4463 

fl AA7A 


0.2297 
0.2401 
0.2417 


Tr 'A t\ ciT nr in a 1 1 on 

procedure 


f40 16 R) 
(60,32,66) 
(120,64,130) 
(240,128,258) 


2.0691 
2.0760 
2.0775 
2.0780 


1.3140 
1.3291 
1.3319 
1.3329 


0.7721 
0.7908 
0.7940 
0.7951 


0.4294 
0.4442 
0.4467 
0.4476 


0.2325 
0.2405 
0.2419 
0.2424 


Our Tree Method 


N 














150 
250 
350 


2.0791 
2.0786 
2.0790 


1.3362 
1.3338 
1.3339 


0.7957 
0.7964 
0.7964 


0.4495 
0.4501 
0.4485 


0.2435 
0.2435 
0.2440 


Reference Value 




2.0784 


1.3336 


0.7960 


0.4483 


0.2428 
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Table 6. Comparison of our method and the semi-closed solution for fixed- 
strike geometric Asian call options for : So = 100, vq = 0.09, r = 0.05, 
k = 1.15, 6 = 0.348, p = -0.64, jj = 0.39. 



T 


K 


MC on Tree with N = 300 


Semi- 
Closed 
Solution 


Difference % 


NumSim 


NumSim 


10 5 


5 * 10 5 


10 6 


10 5 


5 * 10 5 


10 6 


0.2 


90 


10.6598 


10.6551 


10.6562 


10.6571 


0.02 


-0.02 


-0.01 


0.2 


95 


6.6006 


6.5970 


6.5888 


6.5871 


0.20 


0.15 


0.03 


0.2 


100 


3.4699 


3.4564 


3.4510 


3.4478 


0.64 


0.25 


0.09 


0.2 


105 


1.4697 


1.4610 


1.4611 


1.4552 


1.00 


0.40 


0.40 


0.2 


110 


0.4730 


0.4742 


0.4719 


0.4724 


0.14 


0.38 


-0.10 


0.4 


90 


11.7310 


11.7111 


11.7077 


11.7112 


0.17 


0.00 


-0.03 


0.4 


95 


8.0988 


8.1067 


8.0877 


8.0894 


0.12 


0.21 


-0.02 


0.4 


100 


5.1480 


5.1746 


5.1641 


5.1616 


-0.26 


0.25 


0.05 


0.4 


105 


3.0414 


3.0060 


3.0040 


3.0018 


1.32 


0.14 


0.07 


0.4 


110 


1.5555 


1.5776 


1.5679 


1.5715 


-1.02 


0.39 


-0.23 


0.5 


90 


12.2974 


12.2495 


12.2330 


12.2329 


0.53 


0.14 


0.00 


0.5 


95 


8.7711 


8.7668 


8.7753 


8.7553 


0.18 


0.13 


0.23 


0.5 


100 


5.9036 


5.9151 


5.9008 


5.8971 


0.11 


0.31 


0.06 


0.5 


105 


3.7150 


3.7120 


3.7165 


3.7072 


0.21 


0.13 


0.25 


0.5 


110 


2.1622 


2.1692 


2.1595 


2.1589 


0.15 


0.48 


0.03 


1 


90 


14.5646 


14.6087 


14.5937 


14.5779 


-0.09 


0.21 


0.11 


1 


95 


11.6287 


11.5518 


11.5474 


11.5551 


0.64 


-0.03 


-0.07 


1 


100 


8.9708 


8.9378 


8.9530 


8.9457 


0.28 


-0.09 


0.08 


1 


105 


6.8003 


6.7392 


6.7505 


6.7559 


0.66 


-0.25 


-0.08 


1 


110 


5.0161 


4.9878 


4.9704 


4.9722 


0.88 


0.31 


-0.04 


1.5 


90 


16.3889 


16.4588 


16.5200 


16.5030 


-0.69 


-0.27 


0.10 


1.5 


95 


13.7324 


13.7764 


13.7690 


13.7625 


-0.22 


0.10 


0.05 


1.5 


100 


11.3599 


11.3247 


11.3304 


11.3374 


0.20 


-0.11 


-0.06 


1.5 


105 


9.2487 


9.2187 


9.2076 


9.2245 


0.26 


-0.06 


-0.18 


1.5 


110 


7.4342 


7.3959 


7.4019 


7.4122 


0.30 


-0.22 


-0.14 


2 


90 


18.0757 


18.1112 


18.0816 


18.0914 


-0.09 


0.11 


-0.05 


2 


95 


15.6133 


15.6021 


15.5211 


15.5640 


0.32 


0.24 


-0.28 


2 


100 


13.3624 


13.3245 


13.2833 


13.2933 


0.52 


0.24 


-0.08 


2 


105 


11.2855 


11.2862 


11.2627 


11.2728 


0.11 


0.12 


-0.09 


2 


110 


9.4243 


9.4840 


9.4901 


9.4921 


-0.71 


-0.09 


-0.02 


3 


90 


20.6523 


20.4276 


20.5149 


20.5102 


0.69 


-0.40 


0.02 


3 


95 


18.3985 


18.2361 


18.2884 


18.3060 


0.51 


-0.38 


-0.10 


3 


100 


16.2151 


16.2555 


16.2609 


16.2895 


-0.46 


-0.21 


-0.18 


3 


105 


14.5000 


14.4330 


14.4046 


14.4531 


0.32 


-0.14 


-0.34 


3 


110 


12.6065 


12.8177 


12.7982 


12.7882 


-1.42 


0.23 


0.08 
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Table 7. Confidence Intervals for fixed-strike geometric Asian call optk 
for : S = 100, u = 0.09, r = 0.05, k = 1.15, 6 = 0.348, p = -0. 
77 = 0.39. 



Confidence Intervals 95 % 



NumSim = 10 5 


NumSim = 5 * 10 5 


NumSim = 10 6 


(10.6135,10.7060) 


(10.6345, 10.6758) 


(10.6416,10.6708) 


(6.5609,6.6402) 


(6.5793,6.6147) 


(6.5763,6.6014) 


(3.4397,3.5001) 


(3.4429,3.4699) 


(3.4415,3.4605) 


(1.4501,1.4894) 


(1.4522,1.4698) 


(1.4548, 1.4673) 


(0.4623,0.4838) 


(0.4694,0.4790) 


(0.4685,0.4753) 


(11.6678,11.7941) 


(11.6829,11.7394) 


(11.6877,11.7277) 


(8.0438,8.1538) 


(8.0820,8.1313) 


(8.0703,8.1051) 


(5.1027,5.1932) 


(5.1543,5.1948) 


(5.1498,5.1784) 


(3.0065,3.0764) 


(2.9904,3.0216) 


(2.9930,3.0150) 


(1.5308,1.5803) 


(1.5665,1.5887) 


(1.5601,1.5758) 


(12.2270,12.3679) 


(12.2181,12.2808) 


(12.2108,12.2552) 


(8.7094,8.8328) 


(8.7391,8.7944) 


(8.7557,8.7949) 


(5.8516,5.9556) 


(5.8919,5.9384) 


(5.8843,5.9172) 


(3.6735,3.7566) 


(3.6934,3.7306) 


(3.7034,3.7297) 


(2.1305,2.1938) 


(2.1551,2.1833) 


(2.1495,2.1694) 


(14.4642, 14.6650) 


(14.5638,14.6536) 


(14.5619,14.6255) 


(11.5367,11.7208) 


(11.5109,11.5927) 


(11.5186,11.5763) 


(8.8888,9.0528) 


(8.9013,8.9744) 


(8.9272,8.9789) 


(6.7282,6.8724) 


(6.7072,6.7713) 


(6.7278,6.7732) 


(4.9538,5.0784) 


(4.9601,5.0154) 


(4.9508,4.9899) 


(16.2635,16.5144) 


(16.4023,16.5152) 


(16.4800,16.5599) 


(13.6150,13.8498) 


(13.7239,13.8289) 


(13.7319,13.8061) 


(11.2523,11.4676) 


(11.2765,11.3729) 


(11.2963,11.3645) 


(9.1503,9.3471) 


(9.1749,9.2626) 


(9.1766,9.2387) 


(7.3457,7.5226) 


(7.3563,7.4355) 


(7.3739,7.4299) 


(17.9261,18.2253) 


(18.0442,18.1782) 


(18.0342,18.1289) 


(15.4721,15.7544) 


(15.5392,15.6651) 


(15.4767,15.5654) 


(13.2303,13.4945) 


(13.2656,13.3835) 


(13.2416,13.3249) 


(11.1626,11.4084) 


(11.2315,11.3409) 


(11.2240,11.3013) 


(9.3113,9.5373) 


(9.4334,9.5345) 


(9.4544,9.5258) 


(20.4610,20.8436) 


(20.3429,20.5123) 


(20.4547,20.5750) 


(18.2156,18.5814) 


(18.1549,18.3174) 


(18.2310,18.3459) 


(16.0417,16.3885) 


(16.1781,16.3329) 


(16.2062,16.3157) 


(14.3343,14.6656) 


(14.3594, 14.5067) 


(14.3525,14.4566) 


(12.4514,12.7617) 


(12.7476,12.8878) 


(12.7489,12.8476) 
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Table 8. Comparison of our method and the Functional Quantization 
Method by Pages and Printems 25 . Parameters: So = 50, Uq = 0.01, 
r = 0.05, k = 2, 6 = 0.01, p = 0.5, 77 = 0.1. 
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56 


0.125 


(1.04%) 


0.15 


(19.96%) 


0.125 


(2.84%) 


0.1317 


( 0.1303 , 0.1330 ) 
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Table 9. Comparison of our method and Eulcr Simulation for Lookback 
Call option with fixed strike. Parameters: So — 100, v$ — 0.16, r = 0.05, 
k = 3, 9 = 0.04, p = -0.7, r\ = 0.1. 
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time in seconds 




N = 200 


5.71 


Put and Call 


N = 350 


30.37 


(Table 1 & 2) 


N = 500 


89.27 


AniGriccLii Put 
(Table 3) 


N = 250 


13.97 


American 


N = 150 


3.15 


Put 


N = 250 


14.66 


(Table 4 & 5) 


A^ = 350 


40.50 


Geometric 


NumSim = 10 


8.17 


Asian 


NuTTtSivrt — 5 * 10 


40.91 


(Table 6 & 7) 


NumSim = 10 6 


81.79 


Arithmetic 


N = 300 




Asian 
(Tabic 8) 


NumSim = 10 


98.65 


Lookback 
(Table 9) 


N = 3000 

NumSim = 10 


94.54 



